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Abstract. Organisms often grow, migrate and compete in liquid environments, as 
\J^ , well as on solid surfaces. However, relatively little is known about what happens 

Q^ ■ when competing species are mixed and compressed by fluid turbulence. In these 

lectures we review our recent work on population dynamics and population ge- 
netics in compressible velocity fields of one and two dimensions. We discuss why 
1.^ ■ compressible turbulence is relevant for population dynamics in the ocean and we 

consider cases both where the velocity field is turbulent and when it is static. 
Furthermore, we investigate populations in terms of a continuos density field and 
when the populations are treated via discrete particles. In the last case we focus 
on the competition and fixation of one species compared to another 
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1 Introduction 

m_ 

>0 , Challenging problems arise when spatial migrations of species are combined with population 

^T) • genetics. Stochastic number fluctuations are inevitable at a frontier, where the population size 

f^ I is small and the discrete nature of the organisms becomes essential. Depending on the parame- 

O^l ■ ter values, these fluctuations can produce important changes with respect to the deterministic 

predictions [1I2| . When two or more species undergo a Darwinian competition in a spatial en- 
vironment, one must deal with additional issues such as genetic drift (stochastic fluctuations 
in the local fraction of one species compared to another) and Fisher genetic waves, [3] which 
allow more fit species to replace less fit ones. On solid surfaces, the complexities of spatial pop- 
ulation genetics are elegantly accounted for by the stepping stone model, originally introduced 
by Kimura and Weiss |5 j H] ■ 

However, much of population genetics, from the distant past up to the present, played out in 
liquid environments, such as lakes, rivers and oceans. For example, there are fossil evidence for 
oceanic photosynthetic cyanobactcria (likely pre- cursors of chloroplasts in plants and a major 
source of oxygen in the atmosphere) that date back a billion years or nrore [5] . In addition, it 
has recently become possible to perform satellite observations of chlorophyll concentrations to 
identify fluid dynanrical niches of phytoplankton types off the eastern coast of the southern tip 
of South Anrerica |7] , where the domains of the species are largely determined by the tangen- 
tial velocity field obtained from satellite altimetry. In cases such as these, spatial growth and 
evolutionary competition take place in the presence of advecting flows, some of them at high 
Reynolds numbers [5]. 
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Phytoplankton needs light and nutrients to grow and many phytoplankton species are able to 
adjust their density and swim to stay near the surface. Nutrients are brought to the surface from 
deeper ocean layers, usually below 500 meters. Therefore, oceanic circulation plays an important 
role in shaping spatial growth and evolution of plankton species. To appreciate the complexity of 
the problem, it is worthwhile to shortly review our present knowledge on the basic mechanisms 
one should consider. In a three dimensional turbulent flow at high Reynolds number, the velocity 
field is fluctuating over a range of scales [L,r]] where L is the scale of energy pumping in 
the system and 77 = (i^'^/e)^/'* is the Kolmogorov dissipation scale. The velocity field is also 
fluctuating in time. According to Kolmogorov theory, one can define the dissipation time scale 
as T^ = ^Jvje. In the upper oceanic mixed layer , forcing is provided by heat and momentum 
exchange with atmosphere and the observed values [H] of e ranges from Wi^'' cw? j sec" up to 
50cm^/sec"^, which implies 77 e [0.01,2]c7r7 and r^ € [0.01, 300]sec. The phytoplankton size lies 
in the range [10,200]^m with a density difference respect to sea water density in the range 
[0.01,0.1]. Advection of individuals in the ocean should be studied by considering all forces 
acting on them. In particular, because of density mismatch and finite size, individuals are not 
advected as simple Lagrangian tracers [l^ [11] , i.e. the velocity field experienced by each 
individual is not the Lagrangian velocity field, but an effective velocity field which may be not 
incompressible. A suitable measure of compressibility can be defined as k = {(div v)^)/((Vv)^), 
where (..) stands for space and time average. Using the above mentioned values of phytoplankton 
size a, density mismatch ipj p and turbulent energy dissipation e, one obtains 

V^=^^e[i0-,.4] 

P VTr, 

Another very important feature to be considered is the ability of individuals to swim in a pref- 
erential direction towards the largest concentration of nutrients (chemotaxis) . The swimming 
velocity Yc is presently estimated in the range [10, 500]/ir7i/sec. Because of turbulent, individuals 
are subject to external forces which try to change the direction. It is observed that with a char- 
acteristic time B ~ 5sec, individuals try to recover the preferential direction. This mechanism, 
named gyrotaxis [12| and [13) . introduces an effective compressible flow with compressibility 

x/;^ = -5^ e [2.510-3, 1] 
7/ 

It is important to remark that turbulent flows with an effective compressibility can dramati- 
cally change population dynamics: concentration of individuals increases in low pressure regions 
(sinks) and decreases in high pressure regions (source) and the population is spatially charac- 
terized by small scale patchiness. The above discussion shows that intense turbulent activity 
in the oceanic upper layer may introduce non trivial effect, due to compressibility, in the phy- 
toplankton growth and evolution at rather small scale. The same considerations might be 

relevant for large scale motions. Very large scale oceanic circulation (100 — 300^777.) are char- 
acterized by relatively small Rossby number Ro^ defined as Ro = uh /{fL), where uh is the 
characteristic horizontal velocity, order 0.l77T,/sec, / = 10~'*sec~^ is the Coriolis frequency and 
L is the characteristic large scale circulation. For Ro << 1, the velocity field is close to the 
geostrophic balance, meaning that the Coriolis force balances the pressure gradient. Under such 
circumstances, the vertical velocity w is rather small and it can be estimated to be 0.l77T,7?z/sec 
or equivalently few meters/day. The horizontal velocity can be decomposed in the geostrophic 
component Vg and the non geostrophic part Vq where divn'Vg = and divn'Va + 0^11! = 
with divH = dx + dy. According to quasi-geostrophic dynamics, near the surface there exists 
an effective compressible flow acting on time scale order div Vq ^ 10~^sec~^ much longer than 
the longest population growth rate /i ~ 2 10~^sec~"'^. Therefore, at very large scale, population 
dynamics evolves under the advection of an incompressible flow. The above picture changes dra- 
matically if we consider flows at Ro close to 1. Recent numerical simulations as well as direct 
observations [I1],[I1], [IS] have shown that surface density tends to develop sharp horizontal 
gradients (fronts) especially near by the edge of oceanic eddies. Formation of intense fronts. 
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produced by the enhanced filamentation of surface density [T7] [IH], increases the vertical ad- 
vection and destroy geostrophic balance. As a results two important phenomena seem to take 
place in the ocean at relatively large scale (order 10km) [ll]IlD]: regions of relative large and 
positive vertical velocity (upwelling) tends to increase nutrients for phytoplankton providing 
an increase of total biological mass while regions of negative vertical velocity increases concen- 
tration of the phytoplankton population. Frontogenesis, as it is usually named the formation of 
sharp density gradients, can develop vertical velocity up to few millimcter/scc. Consequently, 
the horizontal velocity near frontogenetic regions is characterized by an effective compressibil- 
ity with divH'v ^ 10^^ [TH], i.e. smaller than the population growing rate. The above picture 
suggests the formation of plankton patchincss on scale ranging from IOOtti to 10 — 30km. As 
a tentative conclusion to our short review of phytoplankton in the ocean, albeit the complexity 
of the problem, it seems important to understand the role of turbulent compressible flows in 
population dynamics and population genetics trying to understand, at least in the simplest 
cases, if a new and non trivial phenomenology can be discovered and its relevance to biological 
evolution. 

In these lectures we review our recent work |21 ) .|22 ) .|23j on population dynamics and pop- 
ulation genetics in compressible velocity fields of one and two dimensions, motivated by the 
above discussion. We consider cases both where the velocity field is turbulent and when it is 
static. Furthermore, we investigate populations in terms of a continuos density field and when 
the populations are treated via discrete particles. In the last case we focus on the competition 
and fixation of one species compared to another. 



2 One dimensional case 

In this section we shall discuss some qualitative and quantitative ideas underlying the effect 
of compressible turbulence on population dynamics. We restrict ourself to the one dimensional 
case where most concepts can be discussed using rather simple analytical tools. 
Upon specializing to one dimension, the Fisher equation reads |24j 

dtc + d^iuc)^ Ddlc + ^ic-bc^ (1) 

Equation ([T]) is relevant for the case of compressible flows, where dxU ^ 0, and for the case when 
the field c(x, t) describes the population of inertial particles or biological species. By suitable 
rescaling of c(x.i), we can always set 6 = ^. In the following, unless stated otherwise, we shall 
assume h = ^ whenever /x 7^ and 5 = for /x = 0. 

The Fisher equation for u = has travelling front solutions which can be computed analyt- 
ically: 

cix.t)^ ; (2) 

[1 + Cea;p(-5Mt/6 ± xVmAD/6)]2 

From © we can see that the non linear wave propagates with velocity vp ~ (D/i)^/^ [5], |24j . 
In Fig. ((D) we show a numerical solution of Eq. ^ with D = 0.005, /U = 1 and u = obtained 
by numerical integration on a space domain of size L — 1 with periodic boundary conditions. 
The figure shows the space-time behaviour of c(a;, t) for c{x, t) = 0.1, 0.3.0, 5.0, 7 and 0.9. With 
initial condition c{x,t = 0) nonzero on only a few grid points centered at a; = i/2, c(x,i) 
spreads with a velocity vp ^ 0.07 and, after a time L/vp ^ 4 reaches the boundary. Note that 
the characteristic size of the Fisher'wave interface thichness is order yj D/ fx. 

Let us consider first the case /z = 6 = 0. In this limit, Eq. ([T]) is just the Fokker-Planck 
equation describing the probability distribution P{x, t) = c(x, t) to find a particle in the range 
(x, X -\- dx) at time i, whose dynamics is given by the stochastic differential equation: 

dx I 

— =«(a;,i) + V2^77(t) (3) 

dt 
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Fig. 1. Contour plot of the numerical simulation of eq. [T]with fj, — 1, D = 0.005 and with periodic 
boundary conditions. The initial conditions are c{x, t) = everywhere expect for few grid points near 
L/2 — 0.5 where c — 1. The horizontal axis represents time while the vertical axis is space. 

where r]{t) is a white noise with {'ri{t)r]{t')) = 6{t — t'). Let us assume for the moment that 
u{x,t) = u{x) is time independent and, moreover, let us take u(x,t) = —r(x — xq). Then, the 
stationary solution of (HI is given by 



P{x, t) = A-^exp[-r{x - x^)f/2D] 



(4) 



where A is a normalization constant. P{x^ t) — P{x) is strongly peaked near the points xq 
and ([4]) tells us that P spreads around xq with a characteristic length of order ^ = ^J D/Pq. 
Hereafter, we shall refer to ^ as " quasi-localization length". 

The same argument can be used to study the effect for a more generic turbulent like one 
dimensional field u{x), still time independent. We can identify r' as a typical gradient of the 
turbulent velocity field u. In a turbulent flow, the velocity field is correlated over spatial scale 
of order Vt,/r where v1/2 is the average kinetic energy of the flow. For P to be localized near 
a generic sink at the point xq, despite spatial variation in the turbulent field, we must require 
that the localization length ^ should be smaller than the turbulent correlation scale v^,/ P , i.e. 



DP 



> 2 



(5) 



Condition ([5]) can be easily understood by considering the simple case of a periodic velocity 
field u, i.e. u = Vt,sin{xvt,/P). In this case, condition ([5]) states that D should be small enough 
for the probability P not to spread over all the minima of u. For small D or equivalently for 
large v1/ P , the solution will be localized near the minima of m, at least for the case of a frozen 
turbulent velocity field u{x). 

The above analysis can be extended for velocity field u{x, t) that depend on both space and 
time. The crucial observation is that, close to the sinks Xi of u(x, t), we should have u{xi, t) ^ 0. 
Thus, although u is a time dependent function, sharp peaks in P(x, t) move quite slowly, simply 
because m(x, t) r^ Q near the maximum of P(cc, t). One can consider a Lagrangian path x{t) such 
that a;(0) = xq, where xq is one particular point where u{xq^Q) = and dxu{x,Q)\x=xo < 0. 
From direct numerical simulation of Lagrangian particles in fully developed turbulence, we 
know that the acceleration of Lagrangian particles is a strongly intermittent quantitiy, i.e. it 
is small most of the time with large (intermittent) bursts. Thus, we expect that the localized 
solution of P follows x{t) for quite long times except for intermittent bursts in the turbulent 
fiow. During such bursts, the position where u = changes abruptly, i.e. almost discontinuosly 
from one point, say cc(i), to another point x{t+5t). During the short time interval 5t, P will drift 
and spread, eventually reforming to become localized again near x(t + 5t). The above discussion 
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suggests that the probability P{x, t) will be localized most of the time in the Lagrangian frame, 
except for short time intervals dt during an intermittent burst. 

From ^ we conclude that for large value of D P{x, t) is spread out, while for small enough 
D, P should be a localized or sharply peaked function of x most of the time. An abrupt 
transition, or at least a sharp crossover, from extended to sharply peaked functions J", should 
be observed for decreasing D. 

It is relatively simple to extend the above analysis for a non zero growth rate /x > 0, see 
also |25j for a time independent flow. The requirement ([S|) is now only a necessary condition 
to observe localization in c. For /i > we must also require that the characteristic gradient on 
scale ^ must be larger than /x, i.e. the effect of the small scale turbulent fluctuations should act 
on a time scale smaller than 1/ fJL. We estimate the gradient on scale ^ as ^w(0/C: where 5v{^) 
is the characteristic velocity difference on scale ^. We invoke the Kolmogorov theory, and set 
5v{£) = v^{(^/Lf/^ to obtain: 

In ([HI), we interpret P as the characteristic velocity gradient of the turbulent flow. Note also 
that 5v{^)/^ < P on the average, which leads to the inequality: 

fi<P (7) 

From (O and ([7]) we also find 

a second necessary condition. 

One may wonder whether a non zero growth rate fi can change our previous conclusions 
about the temporal behavior, and in particular about its effect on the dynamics of the La- 
grangian points where u{x,t) = 0. Consider the solution of (HI) at time t, allow for a spatial 
domain of size L, and introduce the average position 

where Z(t) = L dxc{x,t). Upon assuming for simplicity a single localized solution, we can 
think of Xjn just as the position where most of the bacterial concentration c(x, t) is localized. 
We can compute the time derivative Vm{t) = dxm/dt. After a short computation, we obtain: 

v„T,{t) = zl dx{xjn~x)P{x,tf + u{x,t)P{x,t)dx (10) 

where P{x,t) = c{x,t)/Z{t). Note that w,„ is independent of fi. Moreover, when c is localized 
near Xm, both terms on the r.h.s. of (|10p are close to zero. Thus, Vm can be significantly different 
from zero only if c is no longer localized and the first integral on the r.h.s becomes relevant. 
We can now understand the effect of the non linear term in ([1]): when c{x,t) is localized, the 
non linear term does not affect the value of Vm simply because Vm is close to 0. On the other 
hand, when c(x, t) is extended the non linear term drives the system to the state c = 1 which 
is an exact solution in the absence of turbulent convection u{x, t) = 0. 

We now discuss whether our previous analysis can be compared against numerical simu- 
lations of (IT} in the one dimensional case. To completely specify equation (IT} we must define 
the dynamics of the "turbulent" velocity field u(x,t). Although we consider a one dimensional 
case, we want to study the statistical properties of c{x, t) subjected to turbulent fluctuations 
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which are close to those generated by the three dimensional Navier-Stokes equations. Hence, 
the statistical properties oiu(x,t) should be characterized by intermitteney both in space and 
in time. Although intermitteney is not a crucial point in our investigations, we want to use a 
one dimensional velocity field with some generic features in terms of space and time dynamics. 
For this reason, we build the turbulent field u{x, t) by appealing to a simplified shell model 
of fluid turbulence US]- The wavenumber space is divided into shells of scale fc„ = 2'"-^^ko, 
n = 1, 2, .... For each shell with characteristic wavenumber fc„, we describe turbulence by using 
the complex Fourier-like variable Un(t), satisfing the following equation of motion: 

(-77 + '^kl)Un = i(fcn+lU*j+iUn+2 " <^^nU*i_l "n+l 

+ (1 - S)kn-lUn-lUn-2) + In ■ (H) 

The model contains one free parameter, (5, and it conserves two quadratic invariants (when 
the force and the dissipation terms are absent) for all values of 6. The first is the total energy 
Tin l"""!^ ^^^"i t^^*^ second is X)n(~-'^)"^n l'""l^' where a = log2(l— <5). In this note we fix (5 = —0.4. 
For this value of 5 the model reproduces intermitteney features of the real three dimensional 
Navier Stokes equation with surprising good accuracy |26] . Using u„ , we can build the real one 
dimensional velocity field u{x, t) as follows: 

M(x,i) =F^Ke''="^ + <e-"="^], (12) 

n 

where i^ is a free parameter to tunc the strength of velocity fluctuations (given by u„) relative to 
other parameters in the model (see next section) . In all numerical simulations we use a forcing 
function /„ = (e(l -\-^)/u\)5n,l^ i-e. energy is supplied only to the largest scale corresponding to 
n = 1. With this choice, the input power in the shell model is simply given by l/2X]rJ"n/n + 
M„/*] = e , i.e. it is constant in time. To solve Eqs. ([IJ and (jlip we use a finite difference scheme 
with periodic boundary conditions. Theses model equations can be studied in detail without 
major computational efforts. The free parameters of the model are the diffusion constant -D, 
the size of the periodic Id spatial domain L, the growth rate /U, the viscosity v (which fixes the 
Reynolds number Re) , the "strength' of the turbulence F and finally the power input in the 
shell model, namely e. Note that according to the Kolmogorov theory [57], e ^ u^^^/L where 
w^^j, is the mean square velocity. Since Urms ~ F, we obtain that F and e are related as e '--^ F^. 
By rescaling of space, we can always put i = 1. We fix e = 0.04 and ly = 10^®, corresponding 
to an equivalent Re = UrmsL/v ~ 3 x 10^. Most of our numerical results are independent of Re 
when Re is large enough. In the limit Re — ?> oo, the statistical properties of eq. ([T]) depend on the 
remaining free parameters, £*, jj, and F. For future reference, we compare the characteristic time 
scales for this simple model of homogeneous isotropic turbulence with the local doubling times 
of microorganisms described by Eqs. (1). Upon assuming the usual Kolmogorov scaling picture, 
we expect fluid mixing time scales t in the range v^l"^ je^l"^ < t < L^^^/e^^^, or 0.01 <t< 3.0, for 
typical parameter values of the shell model given above. On the other hand, the characteristic 
doubling time ^2 of, say, bacteria, in our model is <2 ~ 1/m- Our simulations typically take 
/x = 1 so that 0.2 < t2 "£ 1.0, implying cell division times somewhere in the middle of the 
Kolmogorov range. Microorganisms that grow rapidly compared to a range of turbulent mixing 
times out to the Kolmogorov outer scale, as is the case here, are crucial to the interesting effects 
we find when /x > 0. Bacteria or yeast, often mechanically shaken at frequencies of order IHz in 
a test tube in standard laboratory protocols, have cell division times of 20-90 minutes, and do 
not satisfy this criterion. However, conditions that approximately match our simulations can 
be found for, say, bacterioplankton in the upper layer of the ocean, where large eddy turnover 
times do exceed microorganism doubling times [28], |29| . 

In agreement with our previous theoretical analysis, in figure ^ we show the numerical 
solutions of Eq. ([1]) for a relatively "strong" turbulent fiow. A striking result is displayed: we 
see no trace of a propagating front: instead, a well-localized pattern of c{x,t) forms and stays 
more or less in a stationary position. For us. Fig. ^ shows a counter intuitive result. One naive 
expectation might be that turbulence enhances mixing. The mixing effect due to turbulence is 
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Fig. 2. Same parameters and initial condition as in Fig. ([1} for equation ([T]) with a "strong turbulent" 
flow u advecting c{x,t). 




Fig. 3. The behavior in time of the total "mass" Z{t) = j dxc{x,t): circles show the function Z for 
the case of Fig. ([T]), i.e. a Fisher wave with no turbulence; triangles show Z for the case of Fig. ([2]) 
when a strong turbulent flows is advecting c{x,t) 



usually parametrized in the literature (27] by assuming an effective (eddy) diffusion coefficient 
Deff ^ D. As a, consequence, one naive guess for Eq. ([T|) is that the spreading of an initial 
population is qualitatively similar to the travelling Fisher wave with a more diffuse interface of 
width ^ De.j f I \i. As wc have seen, this naive prediction is wrong for strong enough turbulence: 
the solution of equation ([T|) shows remarkable localized features which are preserved on time 
scales longer than the characteristic growth time 1/// or even the Fisher wave propagation time 
Ljvp. An important consequence of the localization effect is that the global "mass" (of growing 
microorganisms, say) , Z = J dxc{x, t), behaves differently with and without turbulence. In Fig. 
([3]), we show Z{t): the curve with circles refers to the conditions shown in Fig. ((!])), while the 
curve with triangles to Fig. ([2|). 

The behavior of Z for the Fisher equation without turbulence is a familiar S'-shaped curve 
that reaches the maximum Z = 1 on a time scale L/vp- On the other hand, the effect of 
turbulence (because of localization) on the Fisher equation dynamics reduces significantly Z 
almost by one order of magnitude 

With biological applications in mind, it is important to determine conditions such that the 
spatial distribution of microbial organisms and the carrying capacity of the medium are signifi- 
cantly altered by convective turbulence. Within the framework of the Fisher equation, localiza- 
tion effect has been studied for a constant convection velocity and quenched time-independent 
spatial dependence in the growth rate fi [30], [21], [32], [33] ■ In our case, localization, when it 
happens, is a time-dependent feature and depends on the statistical properties of the compress- 
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ible turbulent flows. It is worth noting that the localized "boom and bust" population cycles 
studied here may significantly effect "gene surfing" [S^ at the edge of a growing population, 
i.e. by changing the probability of gene mutation and fixation in the population. 

One prediction of eq.s ([S]) and ([7]) is that the limit ^ — >■ should be singular. More precisely, 
the quantity {Z{t)) must be equal to 1 for /z = 0, while our predictions based on ([S]) and ([7]) 
imply that {Z{t)) < 1 for /i 7^ because of "quasi localization" of the solutions. In the insert of 
figure (|4|) we show the time averaged {Z{t)), computed for different values of /x for F = 0.5. For 
large fi^ < Z >— >■ 1, as predicted by our phenomcnological approach, while in the limit fi -^ 
the values of < Z > converges to 0.1. To predict the limit /i — > we can assume that c^{x,t) 



by the 
(13) 



for small enough fi can be obtained by the knowledge of the solution co(a;,i) at fi 
relation 

Cf^{x,t) ^^ Zf^co{x,t) 

where in the above equation the symbol = means "in the statistical sense" and Z^ 

(the subscript x indicates average on space). Since the solution Cfj_{x,t) satisfies the constrain 

{ctj.)x - {cl)x = for any fi, we obtain: 



^flf X 



Zlic^ 



01 X 



{cD 



(14) 



Once again we remark that eq. p4|) should be interpreted in a statistical sense, i.e. the time 



average of Z^^ should be equal for small /i to the time average of (cq)^ 
(HI the blue dotted line corresponds to the time average of (cg)^ 



-1. 



In the insert of figure 
: equation ([Tl)) is clearly 



confirmed by our numerical findings. As we shall see in the next section, the same argument 
can be applied for two dimensional compressible flow. 




1*1: 



Fig. 4. Behavior of the carrying capacity {Z)^ as a function of ^r from 128'^ ( dots) and 512'^ (squares) 
numerical simulations with Sc = 1. Note that for ^lt — >■ 0.001, the carrying capacity approaches the 
limit 1/{P'^) (dotted line) predicted by Eq. (|21|) . In the inset we show similar results for one dimensional 
compressible turbulent flows. 



3 Fisher equation in two dimensional compressible flows 



As discussed in the previous section, an advecting compressible turbulent flow leads to highly 
non-trivial dynamics for the Fisher equation. Although previous results were obtained only in 
one dimension using a synthetic advecting flow from a shell mdel of turbulence, two striking 
effects were observed: the concentration field c{x,t) is strongly localized near transient but 
long-lived sinks of the turbulent flows for small enough growth rate fi] in the same limit, the 
space-time average concentration (denoted in the following as carrying capacity) becomes much 
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smaller than its maximuni value 1. Here, we present numerical results aimed at understand- 
ing the behavior of the Fisher'equation for two dimensional compressible turbulent flows and 
extending our previous results to more realistic two dimensional turbulent flows. Our model 
consists in assuming that the microorganism concentration field c{x,t), whose dynamics is 
described by the equation 

dtc + V-{uc) = DV^c + ^lc{l-c) (15) 

We assume that the population is constrained on a planar surface of constant height in a three 
dimensional fully developed turbulent flow with periodic boundary conditions. Such a system 
could be a rough approximation to microorganisms that actively control their bouyancy to 
mantain a fixed depth below the surface of a turbulent fluid. As a consequence of this choice, the 
flow field in the two dimensional slice becomes compressible [35| . We consider here a turbulent 
advecting field u(x, t) described by the Navier-Stokes equations, and nondimensionalize time 
by the Kolmogorov time-scale r = (yje)^^'^ and space by the Kolmogorov length-scale 77 = 
(i/^/e)"'^/^. The non-dimensional numbers charecterizing the evolution of the scalar field C(x, i) 
are then the Schmidt number Sc = v/D and the non-dimensional time fj.T. A particularly 
interesting regime arises when the doubling time fi~^ is somewhere in the middle of the range 
of eddy turnover times that characterize the turbulence. Although the underlying turbulent 
energy cascade is somewhat different [12] , this situation arises for oceanic plankton, who double 
in ^ 12 hours, in a medium with eddy turnover times varying from minutes to months f43| . 

We conducted a three dimensional direct numerical simulation (DNS) of homogeneous, 
isotropic turbulence at two different resolutions (128'^ and 512'^ collocation points) in a cubic 
box of length L = 27r. The Taylor microscale Reynolds number j^T] for the full 3D simulation 
was Rex = 75 and 180, respectively, the viscosities were v = 0.01 and v = 2.05 • 10"'^, the 
total energy dissipation rate was around e ~ 1 in both cases. For the analysis of the Fisher 
equation we focused only on the time evolution of a particular 2D slab taken out of the full 
three dimensional velocity field and evolved a concentration field c{x, t) constrained to lie on this 
plane only. A typical plot of the 2d concentration field, along with the corresponding velocity 
divergence field (taken at time t = 86, Rex ~ 180) in the plane is shown in Fig. [S] {Sc = 5.12): 
the concentration c[x, y, t) is highly peaked in small areas, resembling one dimensional filaments. 
When the microorganisms grow faster than the turnover times of a significant fraction of the 
turbulent eddies, c(x, t) grows in a quasi-static compressible velocity field, and accumulates 
near sinks and along slowly contracting eigendirections, leading to filaments. The geometry of 
the concentration field suggests that c(x,t) is different from zero on a set of fractal dimension 
dp much smaller than 2. A box counting analysis of the fractal dimension of c(x,t) supports 
this view and provides evidence that dp ~ 1. ± 0.15. 

Note that for /i = 0, Eq. (J15p reduces to the Fokker-Planck equation describing the proba- 
bility distribution P(x,y,t) to find a Lagrangian particle subject to a force field u(x,i) at x,y 
at time t: 

dP 

^— + V • (uP) = DV^P (16) 

ot 

The statistical properties of P have been studied in several works (e.g. |33] and [ID]) and it is 
known that for compressible turbulence P(x, i) exhibits a nontrivial multifractal scaling. Upon 
multiplying eqn. ([T6| by P and integrating in space we obtain: \dt {P"^) + h (^^(V • u)) = 
—D ((VP)^) where {■■■)„ denotes a spatial integration. In the statistically stationary regime, 
the above equation reduces to: 

\{P\V-n))=-D{{WPf), (17) 

where now (...) stands for space and time average. Eq. ([T7| shows that for V • u = the 
only possible solution is P = const. However, compressibility leads to nontrivial dynamics 
such that P^ and V • u are anticorrelated. We measure the degree of compressibility by the 
factor K = ((V • u)2)/((Vu)2), and estimate the l.h.s. of Eq. ([TT]) by assuming (p2(V • u)) = 
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— ^i(P^)((V • u)^)^/^, where we used the so called one point closure for turbulent flows p7] 
and Ai is expected to be order unity. We estimate the r.h.s of Eq. (|T7| by assuming: 



{(VPf) = A, 






(18) 



where we define ^ the "quasi-localization" length of P, which is expected to be of the same 
order of the width of the narrow filaments in Fig. [51 Finally we set ((Vu)^) = e/v where e is 
the mean rate of energy dissipation and v is the viscosity. On putting everything together we 
find a localization length given by: 

2A2D^ 



AxJlTe 



e 



One important quantity -from the biological point of view- is the carrying capacity, 



(19) 



Z{t) = 



1 

Z2 



dxdyc{x; t), 



(20) 



and in particular its time average in the statistical steady state with growth rate fj,, {Z(t))^. 
We are interested to understand how {Z)fj_ behaves as a function of ^, in the two important 
limits /x — >■ (X) and ^ — > 0. In the limit // — >■ cxj, we expect the carrying capacity attains 
its maximum value (^)p^oo ~ 1, because when the characteristic time l/n becomes much 
smaller than the Kolmogorov dissipation time r,, = {u/eY^'^, the effect of the velocity field is 
a relatively small perturbation on the rapid growth of the microorganisms. Indeed, consider a 
perturbation expansion of c(x, t) in terms of (5 = 1//.J. On defining c = 2_] <5'ci(x, i), substituting 



in Eq. (|T5|) . assuming steady state, and collecting the terms up to 0{S^) we find c « 1 — e(V • 
u) + e^jV • [m(V • It)] - DW'^{'V • m) - (V • m)^} + 0{d^). The above analysis shows that in 
the limit /i — > oo the concentration field tends to become uniform with the leading correction 
coming from the local compressibility. After substituting the expansion of c in Eq. (j20p one gets 
Z w 1 — (S^/L) J(V • u)'^dx -|- 0{S^). Note that the leading correction to the carrying capacity 
is of order S'^ , is consistent with the physical picture presented above. 




Fig. 5. Plot of the concentration c{x,y,t). The white color indicate regions with low concentration 
while regions of high concentration are denoted by black. 
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By defining F = ((V • tt)^)^/^ as the r.m.s value of the velocity divergence, we expect a 
crossover in the behavior of {Z) ^ for fi < F. In the limit /i — > 0, following our discussion in the 
previous section, we expect that: 

We have tested both Ea. (|2T1) and the limit /i — > oo against our numerical simulations. In Fig. 
(j3]) we show the behavior of (Z)^j for the numerical simulations discussed in this section. The 
horizontal line represents the value 1/(P^) obtained by solving Eq. ([T5| for the same velocity 
field and fi = 0. For our numerical simulations we estimate F = 4.0 and we observe, for ^ > F 
the carrying capacity (^)^ becomes close to its maximum value 1. The limit /z — )■ requires 
some care. Let us define Tb = l//i to be the time scale for the bacteria to grow. The effect 
of turbulence is relevant for Tf, longer than the Kolmogorov dissipation time scale r^ . We also 
expect that Tb must be smaller than the large scale correlation time r^ ~ (L^/e)^/^, which 
depends on the forcing mechanism driving the turbulent flows and the large scale L. Thus, the 
limit /z — >■ can be investigated either for L — >■ cxj or by forcing the system with a constant 
energy input which slows down the large scale, as it is the case in our numerical simulations 

The limit /i — > can be investigated more accurately as follows: according to known results 
on Lagrangian particles in compressible turbulent flows, we know that P should have a multi- 
fractal structure in the inviscid limit i/ —>■ [TT]. If our assumption leading to Eq. (PT|) is correct, 
c(x, t) must show multifractal behavior in the same limit with multifractal exponents similar 
to those of P. For analytical results, see Ref. [331 • Numerical evidence for the multifractal 
behavior of Lagrangian tracers in compressible flows can be found in Ref. [38] . 

We perform a multifractal analysis of the concentration field c{x,y,t) with // > by con- 
sidering the average quantity Cfj_{r) = ^ JBir) ''(^' V^ t)dxdy where B{r) is a square box of size 
r. Then the quantities (^^(r)^) are expected to be scaling functions of r, i.e. {Cf_i{r)P) ^ r'^(p\ 
where a{p) is a non linear function of p with a(2) = —0.47, see [22] for details. 

Our multifractal analysis allow us to investigate the possible relation between the localiza- 
tion length ^ defined in Eq. (|18p and the carrying capacity {Z)^. The localization length ^ can 
be considered as the smallest scale below which one should observe fiuctuations of c{x, t). Thus 
we can expect that (P^(a;,i)) '--^ ^"(2). Using (PT|) we obtain {Z) ^ ^-"(2). In the inset of Figure 
^ we show (Z) as a function of £, (obtained by using ([T5)) for fi = 0.01 and different values of 
the diffusivitiy D. According to Eq. ([19)) . reducing the diffusivity D will shrink the localization 
length ^ and hence (Z)^. From Figure ^ a clear power law behavior is observed with a scaling 
exponent 0.46 very close to the predicted behavior — a(2) = 0.47. 

Finally, we discuss bacterial populations subject to both turbulence and uniform drift be- 
cause of, e.g., sedimentation under the action of gravity field. In this case, we can decompose 
the velocity field into zero mean turbulence fluctuations plus a constant "wind" velocity uq. In 
presence of a mean drift velocity Eg. [T5l becomes: 

(JC 

— + V-[(u + uqcM = D'^'^c + /ic(l - c) (22) 

ot 

where e^ is the unit vector along the x-direction. Note that the mean drift breaks the Galilean 
invariance as the concentration c is advected by the wind, while turbulent fluctuations u remain 
fixed. In Fig. [Sj we show the variation of carrying capacity versus uq for two different values 
of /i and fixed diffusivity D ~ 0.015. We find that for uq < Urms [urms is the root-mean- 
square turbulent velocity) the carrying capacity Z saturates to a value equal to the value of 
Z in absence of uq i.e., quasilocalization by compressible turbulence dominate the dynamics. 
For Mo > Urms the drift velocity delocalizes the bacterial density thereby causing Z — > 1, in 
agreement with the results discussed in [37]. 

4 Discrete population dynamics 

The population dynamics of a single species expanding into new territory was first studied in 
the pioneering works of Fisher, Kolmogorov, Petrovsky and Piscounov (FKPP) |3I24I1| . Later, 
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Fig. 6. Main figure; plot of {Z) as function of a super imposed external velocity ito for /i = 0.01 
(bullets) and \jl = 0.1 (triangles). Inset: log-log plot of (Z) as a function of the localization length ^ 
defined in Eq. H18|l for uq = 0. The slope is consistent with the prediction < Z >~ ^~''"' discussed in 
the text. The numerical simulations are done for /i = 0.01 and diflferent values of D from D — 0.05 to 
D = 0.001. 



Kimura and Weiss studied individual-based counterparts of the FKPP equation [3], revealing 
the important role of number fluctuations. In particular, stochasticity is inevitable at a frontier, 
where the population size is small and the discrete nature of the individuals becomes essential. 
Depending on the parameter values, fluctuations can produce radical changes with respect to 
the deterministic predictions [1I2| . If f{x,t) is the population fraction of, say, a mutant species 
and 1 — f{x, t) that of the wild type, the stochastic FKPP equation reads in one dimension [S]: 



dtf{x,t) = Ddtf + sf{l - J) + JDgfd - m{x,t) 



(23) 



where D is the spatial diffusion constant, Dg is the genetic diffusion constant (inversely propor- 
tional to the local population size), s is the genetic advantage of the mutant and ^ = ^(x, t) is a 
Gaussian noise, delta-correlated in time and space that must be interpreted using Ito calculus 
[S]- In the neutral case (s = 0), number fluctuations induce a striking effect in the population 
dynamics, namely segregation of the two species. One can show that the dynamics of compet- 
ing species in ID can be characterized by the dynamics of boundaries between the / = and 
/ = 1 states of Eq. [23l which perform a random walk. This effect is theoretically predicted by 
Eq. (|23p and confirmed experimentally in the linear inoculation experiments on neutral variants 
of fluorescently labelled bacteria illustrated in Fig. (jS^) [5T]. 

We study the influence of advection on the dynamics of two distinct populations consisting of 
discrete 'particles'. Due to competition and stochasticity, interactions between two populations 
usually drive one of the two populations to extinction. The average time of this event (the 
/ixation time) is a quantity of great biological interest since it determines the amount of genetic 
and ecological diversity that the system can sustain. Studying competition in a hydrodynamics 
context, where both a compressible velocity field and stochasticity due to finite population sizes 
are present, calls for a nontrivial generalization of Eq. (P5|) . One complication is that, because 
of compressibility, the sum of the concentrations of the two species is no longer invariant during 
the dynamics. 

We have overcome these problems through a off-lattice particle model designed to explore 
how compressible velocity fields affect biological competition. Let us consider two different 
organisms, A and B, which advcct and diffuse in space while undergoing duplication (i.e. cell 
division) and density-dependent annihilation (death) , see Fig. [T] Specifically, we implement the 
following stochastic reactions: each particle of species i ~ A, B duplicates with rate ^i and 
annihilates with a rate /x^ni, where n^ is the number of neighboring particles (of both types) in 
an interaction range 6. Let N be the total number of organisms that can be accomodated in the 
unit interval with total density ca+cb = 1. To reduce the number of parameters, we fix (5 = 1/N 



Will be inserted by the editor 



13 



Cell doubling (rate(j): 



Competition/death 
{rate npn): 



Fig. 7. The possible six birth and death processes in the particle model consisting of two species, A 
(red) and B (green). 
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Fig. 8. (a) Experimental range expansion of the two neutra E. coli strains used in Ref. [41 j . but run 
about one day longer (D. Nelson, unpublished). The black bar at the bottom is due a small crack 
in the agar substrate, (b) Space-time plot of the off-lattice particle model with no advecting velocity 
field. A realization characterized by a pattern similar to the experimental one has been selected for 
illustrative purposes, (c) Particle model with a compressible turbulent velocity field. Simulations are 
run until fixation (disappearance of one of the two species); note the reduced carrying capacity and 
the much faster fixation time in (c). Parameters: N = 10^, D — 10~*, ^ = 1. Parameters of the shell 
model are as in 1211. 



as the average particle spacing in the absence of flow. Further, we set p,A = P-b = Ms = fJ-, but 
take HA = /^(l + s) to allow for a selective advantage (faster reproduction rate) of species A. 
We will start by analyzing in depth the neutral case s = and consider the effect of s > in 
the end of the Letter. In one dimension and with these choices of parameters, our macroscopic 
coupled equations for the densities CAix,t) and CB{x,t) of individuals of type A and B in an 
advecting field v{x,t) read 



dtCA=-dx{vcA) + DdlcA+y^CA{l + s~CA-CB) + crAS, 

dtCB^-dxivCB) + DdlcB+fJ-CBil~CA-CB)+(7BC 



(24) 



with (Ta = y/fJ,CA{l + s + CA + cb)/N and (Jb = ^J ^icb{'^ + ca + Cb)/N. ^{x, t) and ^'{x, t) are 
independent delta-correlated noise sources with an Ito-calculus interpretation as in Eq. (j23[) . 

Simulations of the particle model corresponding to (j24p with v = s = Q result in a dynamics 
similar to the one observed in experiments, as shown in Fig.([8jD). In this simple limit, our model 
can be considered as a grand canonical generalization of Eq (|23p . where the total density of 
individuals ca + cb is now allowed to fluctuate around an average value 1 . We fix the following 
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Fig. 9. Average fixation time r/ for neutral competitions in compressible turbulence and sine wave 
advection, as a function of (left) the reduced carrying capacity {Z)f and (right) forcing intensity F 
(small {Z)f in the left panel corresponds to large forcing in the right panel), (left) Red circles and 
blue triangles are particle simulations. Other symbols denote simulations of the continuum equations 
with different resolutions on the unit interval. The black dashed line is the mean field prediction, 
Tf = N{Z)f/2. In (right), only particle simulations are shown and dashed lines are the theoretical 
prediction Tf = tq + c/F based on boundary domains, with fitted parameters tq = 9.5, c = 3.5 in the 
case of the shell model and tq = 16, c = 1.4 in the case of the sine wave. 



parameters: A^ — lO'^, D = 10^'*, /i = 1 and L = 1 where L is a one dimensional domain 
endowed with periodic boundary conditions. With these parameters, the fixation time Tf would 
be ~ 10'* for the one dimensional FKKP equation, and ~ 10'^ for the well-mixed case. 

Introducing a compressible velocity field v{x, t), via the shell model[TT]as shown in Fig.([Bl;), 
leads to radically different dynamics. Individuals tend to concentrate at long-lived sinks in the 
velocity field. Further, extinction is enhanced and the total number of individuals n(t) present 
at time t is on average smaller than N. 

In order to study how a velocity field changes tj, we first analyze two different velocity 
fields: The first is a velocity field v{x, t) generated by a shell model (Eqs. [TT|) of compressible 
turbulence [21] . reproducing the power spectrum of high Reynolds number turbulence with 
forcing intensity F. The second is a static sine wave, v{x) = i^sin(27ra;), representing a simpler 
case in which only one Fourier mode is present, and thus a single sink, in the advecting field. 
In both cases, periodic boundary conditions on the unit interval are implemented. 

Fig.iO shows the average fixation time r/ for s = in the first two cases, while varying 
the intensity F of advection. In the left panel, we plot the fixation times as a function of the 
time-averaged reduced carrying capacity {Z)p, where Z{t) = n(t)/N is the carrying capacity 
reduction, i.e. the ratio between the actual number of particles and the average number of 
particles TV observed in absence of the velocity field. Plotting vs. {Z)f allows comparisons with 
the mean field prediction, Tf = 27V(Z)i?//i, valid for well mixed systems (black dashed line) 
[5]. For the shell model, we include simulations of the macroscopic equations ([M)) with different 
resolutions (256 and 512 lattice sites on the unit interval), obtaining always similar results for 
Tf VS. (Z)f. 

In all cases, the presence of a spatially varying velocity field leads to a dramatic reduction 
of Tf, compared to mean field theory. The fixation time drops abruptly as soon as (Z) < 1, 
even for very small F. 
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